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Abstract. We consider the statistics of the number of nodal domains aka nodal 
counts for eigenfunctions of separable wave equations in arbitrary dimension. We 
give an explicit expression for the limiting distribution of normalised nodal counts 
and analyse some of its universal properties. Our results are illustrated by detailed 
discussion of simple examples and numerical nodal count distributions. 
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1. Introduction 

We consider real square-integrable eigenfunctions $(q) of the stationary Schrodinger 
equation 

$(q) = -A M $(q) + V(q)$(q) = £$(q) (1) 

for a massive point particle on an s-dimensional smooth connected Riemannian manifold 
A4 with local coordinates q = (q 1 , . . . , q s ) . Here, Am is the Laplace-Beltrami operator 
on A4, V(q) is a potential, and is an energy eigenvalue. We have set the value of the 
physical constant of Planck's constant squared over twice the mass of the particle 
equal to one by appropriate choice of units. 

We will allow that M. has a boundary and will impose boundary conditions on $(q) 
such that the Schrodinger operator H defined in (1) is self-adjoint (e.g. Dirichlet or 
Neumann boundary conditions). 

We consider only non-negative potentials for which the classically allowed region 
Ve — {q : y(q) < E} is compact and connected. This ensures a discrete and non- 
negative energy spectrum. If V(q) = (free motion) the condition implies that the 
manifold M. is compact. 

We arrange the spectrum in ascending order as < E\ < Ei < ■ ■ ■ < En < 
En+i < • • • and denote by <I?jv(q) the eigenfunction corresponding to En. For a given 
eigenfunction $Ar(q) the nodal set Af = $^ 1 (0) C M consists of all points on the 
manifold where the eigenfunction vanishes. A nodal domain V C M of $jv(q) is a 
maximally connected region where the sign of $7v(q) does not change. 

The characterisation of eigenfunctions in terms of their nodal set has a history 
which is more than 200 years old with the first systematic treatment by Chladni [1] 
who visualised the vibration modes of plates with sand that accumulates at the nodal 
set. Among other things he also counted the number of different nodal domains for each 
vibration mode and used these number to characterise the modes for a given shape. The 
number of nodal domains of eigenfunctions or nodal counts will also be the subject of 
the present contribution. For the wave function $ w we denote the nodal count by v N . 
The collection of all nodal counts forms the nodal sequence {z^v}^ =1 (for systems with 
degenerate eigenvalues this definition of the nodal sequence is incomplete). 
In one dimension Sturm's oscillation theorem [2] states vn = N under very general 
conditions. The generalisation of this seminal result to quasi one dimensional systems 
such as quantum graphs has been a recent research topic [3, 4]. For arbitrary dimension 
a seminal result is Courant's nodal domain theorem [5] which states un < N for the 
Laplacian in any dimension. Pleijel later showed that in dimension d = 2 the upper 
bound v N = N is achieved only a finite number of times [6]. 

In recent years it has been established that the nodal sequence contains a lot of 
information about the underlying geometry. It has been conjectured in [7] that the 
nodal count sequence in some cases allows a full reconstruction of the manifold M. up 
to an overall scaling factor, and that it can be used to distinguish between isospectral 
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partners. These conjectures have been partly confirmed and refined in recent years 
[8, 9, 10, 11, 12, 13]. 

Another recent line of research focusses on the statistics of the nodal counts. To 
this end one defines [14] the normalised nodal count by the ratio 

and focusses on the distribution of its values in an energy window. Courant's theorem 
implies < £jv < 1. For a given spectral interval I g (E) = [E, (1 + g)E\ (where g > 0) 
one defines the nodal count distribution formally by 

PWi^J^- E ^-^) ( 3 ) 

^> N: E N e I g (E) 

where Nj i E } is the number of eigenvalues in I g (E). An interesting question concerns 
the existence and properties of a (smooth) limiting distribution 

P(0 = lim P(0i g{E) • (4) 

hi— >oo 

Such a limit may exist (in the weak sense) because the number of states in the interval 
I g {E) grows as E — > oo. 

For two-dimensional separable systems a semiclassical theory shows [14, 15] that the 
limiting function indeed exists and that it can be expressed explicitly in terms of the 
corresponding integrable classical dynamics. It has a number of universal features: 

(i) The limiting distribution has support < £ < £ cr i t where the critical value is smaller 
than one (which is consistent with Pleijel's theorem [6]). 

(ii) Near the critical value the limiting distribution has a square-root singularity 

P(0~(£crit-0" V2 for£<£cnf (5) 

In this work we will generalize this theory to separable systems in any dimension. 

For non-separable systems the semiclassical theory breaks down - mainly due to 
the lack of an explicit functional that maps a given eigenfunction to its nodal count. In 
this case one may still find the nodal count numerically using for instance variants of the 
Hoshen-Kopelman algorithm [16]. For two-dimensional systems with a corresponding 
classical dynamics that shows chaos (this is usually referred to as quantum or wave 
chaos) such an approach revealed that the limiting distribution is universal. Independent 
of the details of the system the limiting distribution contracts to a Gaussian located 
at a universal value £ u , i.e. = 5(£ — £ u ) [14]. Consistency with Berry's random 

wave conjecture [17] has also been checked numerically - the conjecture states that 
eigenfunctions of a chaotic billiard follow the same statistics as the (monochromatic) 
Gaussian random wave model (a random superposition of plane waves of the same 
wavelength) . 

The universality of the nodal count statistics for wave-chaotic systems in two dimensions 
has been explained in a seminal work by Bogomolny and Schmit [18] who constructed 
a heuristic parameter-free critical percolation model that predicts the numerical value 
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of £ u in perfect agreement with numerical calculations and with the Gaussian random 
wave model (see also [19, 20, 21]). Proving rigorously the implied conjecture that the 
two-dimensional Gaussian random wave model and wave functions of chaotic billiards 
are realisations of critical percolation is certainly one of the most challenging open 
mathematical questions in the field. Indeed a few of the implied properties have already 
been proven for random spherical harmonics [22]. A related and equally challenging 
conjecture states that the nodal lines for such systems are a realisation of stochastic 
Loewner evolution (SLE) [19, 21, 23, 24, 25]. The theoretically known statistical 
properties of nodal counts in two-dimensional wave-chaotic systems have also been tested 
thoroughly in experimental settings [26, 27, 28]. 

Preliminary theoretical and numerical results for two-dimensional systems that are 
neither separable nor nor fully wave chaotic have been obtained for non-integrable 
systems with mixed phase space [29] and for integrable systems for which the wave 
equation is not separable [30] . Especially the latter shows that nodal count statistics in 
non-separable integrable systems have a high degree of complexity with a few features 
that resemble either the separable or the wave-chaotic case while new features appear. 

In this work we address nodal counts in arbitrary dimensions. Indeed little is 
known for dimension larger than two. We will focus on the separable case. In Section 
2 we derive an asymptotic expression for the normalised nodal count and related it 
to the geometry of the unit energy shell in action space. In Section 3 we will give a 
general expression for the limiting nodal count distribution and show that it has some 
universal properties whose details change with the dimension. In Section 4 the cuboid 
and the harmonic oscillator are discussed in more detail and the limiting distribution is 
compared to numerically obtained histograms for finite energy intervals. Eventually we 
will discuss in Section 5 some generalisations of our results and also comment on the 
nodal count for wave chaotic systems and random waves in higher dimensions. 

2. Nodal domain distributions for separable systems 

We consider nodal counts for solutions of the wave equation (1) in the case where a 
separation Ansatz leads to the full solution of the eigenvalue problem. The tools we 
will apply for the derivation of the nodal counts and of the limiting distribution (4) are 
EBK quantisation and Poisson summation. The asymptotic limit will be an integral 
over a region in phase space, and it will involve only classical quantities. We will start 
with introducing the relevant classical mechanics. We will not try to be as general as 
possible during the derivation. Rather we will make some assumptions that simplify the 
derivation and later discuss (see Section 5) which assumptions are essential and which 
may be relaxed. 



On the Nodal Count Statistics for Separable Systems in any Dimension 5 
2.1. EBK quantisation and its implication for nodal counts 

Separability of the wave equation (1) implies that there exist coordinates q = (q 1 , . . . , q s ) 
which (almost) cover the whole s-dimensional manifold M, such that any eigenfunction 
can be written in a product form 

s 

*(q) = n* (0 («) • ( 6 ) 

1=1 

For such systems semiclassical Einstein-Brillouin-Keller (EBK) quantisation can be 
performed successfully. The corresponding classical Hamiltonian Mechanics on the phase 
space T*M (cotangent bundle to the configuration manifold M) is generated by the 

s 

#(p,q)= Y.9 UV (<tiPuPv + V{q) (7) 

u,v=l 

where pi is the conjugate momentum to q l , and g uv is the inverse to the metric tensor 
g uv which defines the squared distance ds 2 = v=1 g uv (q)dq u dq v . 
Quantum separability implies that the corresponding Hamiltonian dynamics is 
integrable. The dynamics is confined to an s-dimensional sub-manifold defined by s 
independent constants of motion C n (p, q) = c n in phase space that (generically) has the 
topology of a torus. One introduces the action variables 




where the integration is a long the curve in the pi-q plane where it intersects with the 
torus defined by the values c for the constants of motion - the action is proportional 
to the area enclosed by the torus in that plane. One may perform a canonical 
transformation to action and angle variables (p, q) i— > (1,0) where = (9 l ,...,9 s ) 
are conjugate to the actions I = (I\, . . . ,I S ). I.e. the phase space is foliated in tori 
such that a point in phase space is specified by the torus with action variables I and 
the position on the torus specified in terms of the s angles 0. As the action variables 
are constants of motion all angle variable become cyclic variables for the transformed 
Hamilton function Hil). 

We will make the following additional assumptions on the classical Hamiltonian 
dynamics: 

(Al) The potential is non-negative and the classically allowed region Ve = 
{q G M : V(q) < E} is connected and compact. We have stated this assumption 
in the introduction. This assumption ensures we have a discrete non-negative 
spectrum. 

(A2) There is a one-to-one correspondence between tori in phase space and points 
I in action space. This assumption excludes double-well potentials and similar 
potentials in higher dimensions where action variables can only be defined locally 
in regions bounded by stationary points and separatrices. 
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(A3) A related assumption is the Hamiltonian is a strictly increasing function of all action 
variables 

u\l) ee °— > . (9) 

Here is the angular velocity of the angle variable 9 l on the torus defined by I. 
We also assume that the Hessian matrix j^§fr at any point is non-negative. 

(A4) Each action takes positive values l\ > and is not bounded from above. This 
assumption excludes that any of the variables q l in which the wave function 
separates is cyclic. This is less restrictive than it may appear: for a system with 
rotational invariance one may reduce the attention either to functions which are 
even or odd under a reflection with respect to a hyperplane through the axis of 
rotation. 

(A5) We assume that the Hamilton function is a homogeneous function of the actions. 
For A > we then have 

H(Xl) = X a H{I) (10) 

where a > is the degree of homogeneity. This assumption implies that the 
dynamics in each energy shell is equivalent up to a scaling factor. For free motion 
on a manifold one has a = 2, so this assumption is mainly a restriction on the 
potentials. Note that the harmonic oscillator in any dimension has degree a = 1. 

The above assumptions are not completely independent. Some may be relaxed without 

distorting our discussion too much (see Section 5). 

The EBK spectrum of semiclassical energy eigenvalues is given by 

£™ K = H(h =n 1 + n 1 ,...,I a = n, + n.) (11) 

where the quantum numbers n\ = 0, 1, 2, . . . are non-negative integers, and the shifts 
Hi are fixed numbers of order unity. E.g. the s-dimensional harmonic oscillator has 
fii = 1/2 for all / and free motion on an s-dimensional cuboid with Dirichlet boundary 
conditions has //; = 1. For our discussion the actual value of /i; is not relevant. 
For a given set of quantum numbers the wave function can be written as 

$n(q) = n^ ) (^) (12) 

1=1 

with real functions 4>n\q l ) of one variable. The corresponding nodal pattern will then 
have a checker board structure. Each of these functions obeys Sturm's oscillation 
theorem, i.e. <pn\q l ) contains n\ nodal points. For the explicit EBK wave functions 
this is straight forward to show. This implies that the number of nodal domains in the 
wave function $ n (q) is equal to 

s 

v n = ]J(ni + 1) . (13) 
i=i 

Note that for a degenerate spectrum separability implies a definite and natural choice 
of preferred basis functions. 
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2.2. The normalised nodal counts and Weyl's law 

In order to find the normalised nodal count £jy = Vn/N for a given wave function 
with quantum numbers n = (ni, . . . ,n s ) we need to know the spectral counting index 
N = N(\). An exact ordering of the quantum numbers is a formidable task - in the 
degenerate case one also needs to make some choice for the order of basis functions 
with the same energy. In the present context any such order would be fine - as it 
turns out to leading order one only needs a sufficiently good approximation to the exact 
counting index as provided by Weyl's law. Indeed the semiclassical approximation 
we use introduces an error in the ordering which may easily exceed any influence of 
degeneracies. Weyl's law states that 

N(E) ~ V r E s/a (14) 

gives the leading asymptotic order of the spectral counting index as E — > oo. Here 



Vi 



= J dl (15) 



is the volume of the region T = {I : < H(I) < 1} in action space. For a free particle 
it is related to the volume Vm of the manifold by Vr = Vm where Vb s = r (C^ i s 
the volume of the s-dimensional unit ball. 

Weyl's law (14) allows us to write the asymptotic expression 

e n ~ nr = ifo + i) „ nk^i + 0{E -i /a) (16) 

for the normalised nodal count. The error estimate on the right side of (16) is based 
on the homogeneity of the Hamilton function which implies ni ~ E x l a . Expression (16) 
will serve as the starting point of the derivation of the limiting distribution in section 3. 

Let us now derive Weyl's law (14) in the present setting. This will not only serve 
us to estimate next-to leading orders in the asymptotic formulas (14) and (16) but also 
give us an opportunity to introduce further details of the setting. The exact spectral 
counting function is defined by 

N(E) = #{£ n < E} = ©(£ " E n ) (17) 

n 

where Q(x) is Heaviside's unit step function. Replacing the exact energies by the EBK 
approximation introduces a small error by shifting the positions of the steps slightly. 
The error introduced by this shift is much smaller than the fluctuations in the spectral 
counting function around its mean value and will be neglected. The Poisson summation 
formula in the form 

J2F(n)= / e 2 * Mx F(x)dx (18) 

and the homogeneity of the Hamilton function asserted by assumption (A5) allow us to 
write 



N(E) = E S ' a I e 2^V«M.I-2 7 riM^ e(l _ H ^ d sj 



(19) 
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N(E) + N osc (E) (20) 

where 

N(E) = E s/a (^j d s I + 0(E- 1/a )^j ~ E s/a V r (21) 

is the contribution from the non-oscillating integral Mi = M 2 = . . . = M s = and 
N OS c(E) is the sum over all remaining (oscillating) integrals - each being at most of 
order _E( s ~ 1 )/ a . Altogether we have derived Weyl's law (14) and estimated that the 
sub-leading correction is a factor of order E l l a smaller than the leading term. 



2.3. The geometry of the energy shell and rescaled actions 

It is worth looking at the geometry of the region T and the hyper-surface dF in more 
detail (see Figure 1 for an illustration). Indeed we here deal integrals of the type 



1= / /(IK/ (22) 
Jet 

with a homogeneous function /(I) (of order (3). Here E ■ F — {I : H(T) < E} is 
a scaled version of the region F. Note that F is compact and convex - compactness 
follows from assumption (Al) and convexity from the second part of assumption (A3). 
Indeed, compactness of the classically allowed region Ve implies compactness of the 
region W = {(p, q) : H(p, q) < 1} in phase space because the allowed momenta for any 
point q G V(E) form a closed s-dim ball in the cotangent space T*.M q . Describing the 
region W in action-angle variables eventually implies compactness. 
Using homogeneity one may reduce the s-dimensional integral (22) to an s — 1- 
dimensional integral over the s — 1-dimensional compact surface (unit energy shell) 
dF = {I : H(l) = 1} in momentum space. Note that dF is the non-trivial part of 
the boundary of the region F and intersects the hyperplanes l\ = (which are also 
boundaries of F). 

This reduction is performed by a substitution to rescaled action variables. The latter 
are defined by 

I = eJ (23) 

such that 

H(3) = 1. (24) 

The rescaled action variables J\ are not independent. Assumption (A3) allows us to use 
the implicit function theorem and solve (24) for 

I s = Z r (Ji, . . . , J 9 -i) ■ (25) 

We will denote the s — 1-tuple of rescaled actions that appears as the argument by 
= (J 1; . . . , Jg_i) such that J = (J^, Z r (Jn)). Again using assumption (A2) one can 
show that the function Zr(Jn) is a decreasing of all arguments because 
dZ v ui(J a ,Z r (J n )) 



dJi u s (Jn,Z r (3 n )) 



(26) 
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Figure 1. Illustration of the region Y = {I : H(l) < 1} in s-dimensional action space. 
The s — 1-dimensional hyperplane I s = is represented by a two-dimensional plane 
spanned by the I\ and I s -i axes in this picture. The illustration also shows the hyper- 
surfaces dY (the unit energy shell, i.e. the level set of the Hamilton function H(T) = 1) 
and H These are the upper and lower parts of the boundary of Y. 
A general point in action space with coordinates I is projected onto the unit energy 
shell dY where it is represented by the s-tuple J. 

by the implicit function theorem. The intersection of dT with J s = Zp(Jr>) = marks 
the boundary of the range Q of the variables J^- 
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Let us now come back to the transformation (23). It implies a change of integration 
variables from the s unsealed actions I to the independent s — 1 scaled actions with 
values in the region Q and a scaling factor e G [0, E\. The Jacobean can be calculated 
straight forwardly and is given by 

J = 6 s - 1 (Zr(J n ) - Jn • V Js2 Z r (J n )) = s s - 1 —- (27) 

U) S {J ) 

where the right hand-side follows from Euler's homogeneous function theorem for the 
Hamilton function. The Jacobean is thus positive for Jq E and e E (0,Ej. With the 
shorthand 

s-l 

dV = (Z r (Jn) - Jn • V Jn Z r (J n )) J] rfJ, (28) 

i=i 

we may now rewrite (22) as 

x= e^_ r rfr/(J ^ z(Jn)) (29) 



For / = 1 this implies 



Vr = - / dT- (30) 



s 

It is worth giving a geometrical illustration of the asymptotic nodal count (16). In 
rescaled action variables n\ = 1\ — eJ\ the normalised nodal counts becomes a ratio 



~ (31) 



where 



V(J n ) = \ I[ J iJ Z ( J n) (32) 

is the volume of an s-dimensional cuboid in action space with faces parallel to the 
hyperplanes l\ = 0, and with one corner in the origin and the other on a point 
J = (Jq, Z(Jq)) on the surface dT (see figure 2 for an illustration). As V(Jq) < Vr 
we immediately obtain £(Jq) < 1 which is consistent with Courant's theorem [5]. Since 
the maximal value of the volume V(Jn) is definitely smaller than Vr the result is also 
consistent with Pleijel's theorem [6]. Let J crit be the values for the rescaled action where 
V(j7n) takes its maximal value for Jq G Q. Then J crit is a solution of the equations 

Z(Jn) = -J,^^ Z = l,...,a-1. (33) 

OJi 

Note that the left hand side is a strictly decreasing function of J t while convexity of T 
implies that the right hand side is an increasing function. As a consequence the solution 
to equation (33) is unique and V(J^) has only one critical point in fl which is the global 
maximum. 

In the asymptotic regime E — >■ oo there will be no normalised nodal counts which are 
larger than the critical value 

Cent = < 1 . (34) 
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i 



or 




h 



Figure 2. Illustration of the geometric interpretation of the normalised nodal count 
as a ratio of two volumes £(J$i) = V(Jjj)/Vr < 1- The cuboid of volume V(JnJ is 
inscribed in the region T with volume Vr- The faces of the cuboid are parallel to the 
hyper-surfaces = with one corner at the origin and the opposite corner on the 
surface dT. 

Simple geometric intuition based on this picture shows that £ cr i t will usually not be very 
close to either zero or unity for moderate dimensions - in high dimensions one may have 

£crit < 1. 
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3. The nodal count distribution and its universal properties 

Let us now consider the nodal domain distribution (3). Poisson summation (18) and an 
application of (29) then gives P(£)i g (E) ~ P(i) + 0(E~ l / a ) with the limiting distribution 



The limiting distribution above does not depend on the size g of the spectral interval 
I g {E) = [E, (1 +g)E]. Note also that P(£) is obtained as a weak limit, i.e. in the sense 
of weak convergence of linear functionals which (together with the fact that the support 
is always finite) ensures convergence of all moments. 

In practice one may consider P(£)i (e) numerically in the form of a histogram (i.e. in a 
locally averaged form) and these will have some corrections to the limiting distribution 
- these corrections will depend on the energy E, the size g of the spectral interval, and 
the bin size that has been used for the histogram. As E — > oo with g and bin size fixed 
the fluctuation will become smaller. Indeed one may decrease the bin size moderately 
as E increases - for convergence to a smooth function one just has to ensure that the 
number of normalised nodal counts per bin increases indefinitely. 

Expression (35) is quite general and we will now turn deriving some universal 
properties by a close analysis of this expression. In section 2.3 we have already mentioned 
that there is an upper bound £ crit < 1 to the normalised nodal count. This implies a 
cut-off for the nodal domain distribution P(£) which has its support inside < £ < £ cr it. 
Within its support P(£) is differentiable. This follows from the fact that V(Jn) has only 
one critical point (maximum) in Q. At £ = and £ = £ cr it the distribution P(£) may 
have singularities. We will show that the behaviour at the cut-off is mainly governed by 
the dimension s. For s = 2 one has a square root divergence, for s = 3 there is a finite 
step, and for s > 4 the distribution becomes continuous at £ = £ cr it but not smooth. 

3.1. The behaviour o/P(£) near the cut-off £ cr i t . 

For £ smaller and close to £ cr i t the contribution to P(£) depend on the behaviour of 
V(Jn) near its maximal value which it takes at J cr i t . Taylor expansion of V(Jn) to 
second order around the maximum gives 




(35) 




s-1 



(36) 



where AJi = J\ — J CT it,i and %ui is a positive definite matrix. From 



s 




(37) 



i=i 



one obtains 



p(0 



Z(Jcrit)V l5s -2 




2Vrv / det 7 W 
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where Vss-2 = 2n^ s ~ 1 ^ 2 /T((s — l)/2) is the volume of the s — 2-dimensional sphere. 
The two dimensional case s = 2 is included in this analysis, with V50 = 2. In this 
case P(£) diverges oc /e 1 _ f as was shown before in [14]. When s = 3 we observe that 

V Scrit S 

P(£) ->■ const > as £ ->■ £ crit from below. For s > 4 we have P(£) oc (£ C rit-0 (s ~ 3)/2 ~> 
such that P(£) is continuous at £ = £ crit . 

3.2. The behaviour of P(£) near £ = 0. 

Now, we will study the behaviour of P(£) near £ = 0. For s = 2 it is not difficult to 
show that 

1, ZUA-hZ'UA 1, MZ)-ZJ'AZ) , , 

& P «) = 2 jB, Z( J.) + J, g ( J,) + 2 & J,(z) + z JKZ) = 1 (39) 
where Ji(Z) is the inverse function of Z{J\). 

For the rest of this section we keep our focus on s > 3. Note that the 5-function 
$(£ — V(Jn)/Vr) in expression (35) for £ < £ crit reduces the integral to an s — 2- 
dimensional integral over the level surfaces of V(Jsi). These are closed deformations 
of an s — 2-dimensional sphere. For our present purpose it is useful to write 

dT = sdT 1 + oT 2 (40 a) 

s-1 

dT 1 = Z{3 a )J{dJ l (406) 
1=1 

dT * = ' §TT ^ • V ^ V( Jn)] n dJ, (40c) 

We will show below that dT 2 does not give a contribution to the nodal count distribution 
which then reduces to 

where dS^ is the surface volume (area) element of the surface 

S e = {Jn : V(Jn) = £V r }. (42) 
In order to show that the corresponding integral over aT 2 vanishes one may start with 

where n = Vj n V(Jn)/|Vj n V(J^)| is the unit normal vector on the surface S^. GauB' 
theorem turns this into a volume integral over the region V(Jn) > £Vr enclosed by 
the surface. The corresponding integrand is the divergence of the vector J " which 
vanishes identically which proves that equation (41) is correct. 

For £ — > + the surface S% approaches the boundary dQ of Q. From (41) we see 
that the contributions from a volume element dS$ carry a weight Z(Jn)/|Vj n V(Jn)|, so 
it will be dominated by any critical points where |Vj n V(Jn)| = close to S^. Indeed 
there are such critical points on the boundary dQ and they coincide with the set of 
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points where dQ is not smooth. These are the s — 3-dimensional intersections of any 2 
hyperplanes Ji = or of one such hyperplane with Z(Jq) = 0. All of these are saddle 
points. 

For s = 3 the saddles are isolated at the three corners of dfl. The Hessian at these 
saddle points is not degenerate. 

For s > 4 the saddles are not isolated and the Hessian is degenerate. The saddles form 
continuous surfaces and where they intersect the suppression |Vj n V(Jn)| close to the 
intersection is enhanced by the combined effect of two or more intersecting saddle point 
surfaces. 

The s corners J*^ (c = 0, . . . , s — 1) of f2 thus dominate the P(£) for £ — > + (see figure 
3 for an illustration). Explicitly the corners are given by the origin and the s — 1 
points j( c ) where Z(3q) = intersects with the s — 2 hyperplanes of the form J; = 
where / G {1, . . . , s — 1} — {c}. At the corners the lowest order of a derivative which 
does not vanish identically is (s — 1). We will show that this leads to a divergence of 
the nodal count distribution lim^ + P(0 — °°- 

We may focus on the leading contribution from the corner at the origin which dominates 
the distribution for small £ - indeed the weight Z(J^)/|Vj n V(Jn)| suppresses the 
contribution at the corners due to the factor Z(Jn) which is zero for all corners apart 
from the origin. 

In order to derive the contribution from the origin let us start by expanding enumerator 
and denominator of the weight Z(J^)/|Vj S2 V(Jn)| independently. For the denominator 
one has |Vj n V(J^)| = 0(AJ^ 1 ). The enumerator Z(Jn) however remains finite 
Z(Jq) = Zq > near the origin. Now consider the contribution 

p( ^t//f ? -^n ij ' (44) 



i=i 



from a small region C that contains the origin. We have used Z(Jq) ~ Z with 
corrections O(Jq). The calculation is simplified if we take C as an s — 1-dimensional 
cuboid with side lengths ai 

C = {J n :0< J,<o,,Z = l,...,s-l}. (45) 

The actual values of the side lengths ai will not enter the leading order which implies 
that we have a true corner phenomenon and that the leading asymptotic order of integral 
does not depend on the details of region C. 
In order to perform the integration set 

C = (46) 

and transform coordinates Jq — > (Ji, ■ ■ ■ , J s -2, 0- The ^-integral can be performed and 
leaves 

J i Jo J s-2 \ a s -iZ [[ l=1 Ji ) 
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Figure 3. Illustration of the graph of the function V(r) over the region fl. For s = 3 
this illustration is exact and one can see three saddles in the corners of the region 
il. Near these corners, especially near the one at the origin, the values of V(T) are 
strongly suppressed. 

For s > 3 one the illustration has to be taken with some care as it represents an s — 1- 
dimensional plane by a two-dimensional - the strong suppression is actually enhanced 
in this case. 

where the factors J^ 1 stem from the Jacobean. This integral can be solved iteratively 
using 

a d%n w^/ ,x ^/ 1N (loga) m - (log6)' +1 . . 

—(logx) l Q(x-b) = G(a-b y 6 ' — — (48) 
o x I + 1 
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and gives the leading contribution 

p «> ~ h^r ' (49) 

Note that the leading order corrections to this depend on the lengths a\ which implies 
that a global approach is necessary to evaluate the next-to-leading order of P(£) as 
£ — > + . We see that the leading order diverges as an s — 2-th power of a logarithm 
with a universal constant l/(s — 2)!. Any system dependent features can only enter at 
next-to-leading order. 



3.3. Monotonicity 

For s = 2 one can check directly that P(£) is a strictly increasing function for 
< £ < £ cr i t . Indeed (41) is valid for s = 2 and gives 

prn = gCV) z(j lt+ ) 

lU Z(J 1) _) + Z'(J 1 ,_)J li _ Z(J 1>+ ) + Z'(J 1)+ )J 1>+ 1 ) 

where J\ _ < Ji i+ are the two solutions of £ = JiZ(Ji). Note that both terms in (50) 
are positive. Differentiation of the first term yields 

Ji,-Z'(Ji,-) 2 - Z(Jx,_)Z'( J lf _) - ^(J^Z'Vi,-) d^i,- 



> (51) 

because > 0. Analogously the derivative of the second term gives a positive 

contribution because ^± < 0. 

For s > 3 our calculations above imply that P(£) is a decreasing function in a 
neighbourhood of £ = 0. For s > 4 we have also shown that P(£) is a decreasing 
function near £ = £ crit (for s = 3 our results are consistent with a decreasing function). 
This suggests that P(£) may be a decreasing function over its full support < £ < £ cr it 
for s > 3. Such a conjecture is supported by all example calculations that we have 
performed - however we have not been able to prove it. 

4. Two simple examples: the harmonic oscillator and the cuboid 

Let us now illustrate our results with a few examples that allow for more explicit 
treatment. 

4-1. The s- dimensional harmonic oscillator 

For a harmonic oscillator the Hamilton function is linear in the action variables 

s 

H(I) = J2^i (52) 
i=i 

The unit energy shell dT in action space is then described by the function 

1 5-1 

J s = Z(3 n ) = l YW,. (53) 
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The volume of the region T is Vr = s! ^ ^ • 

From V(Jn) = Z(J n ) YliZi Ji one finds its maximum value at J CT it,i = ^ such that 
V(Jcrit) = jrj=p — — ■ This implies the critical value 



£crit,s = ~ s (54) 

for the normalised nodal count. Note that the individual frequencies do not enter. 
Indeed the complete nodal count distribution with s degrees of freedom does not depend 
on the frequencies and can be expressed as 

s-l 

P s (0 = s\l \l-yj l \x (55) 



s-l \ s-l 



Juzi h<i V i=i , 

x s - S \ (i - j2 j j n j ij n dji ■ ( 5e ) 

For s = 2 this integral has the explicit form 

P 2 (£) = (l-20~ 1/2 for e< 1/2. (57) 
For arbitrary s one may evaluate all positive integer moments 

See Figure 4 for the graph of the limiting distribution (56) for s = 2, 3,4 together with 
numerical data obtained for finite energy intervals. 

4-2. The s- dimensional cuboid 

We consider the free particle in an s-dimensional cuboid (rectangular box) with side 
lengths a\ (I = 1, . . . , s). With Dirichlet conditions on the boundary of the box one 
obtains exact energy eigenvalues 

E » = - 2S L4 (59) 
1=1 1 

where the quantum numbers n\ run over positive integers. The corresponding classical 
Hamilton function H(T) = tx 1 Xw=i is homogeneous of order a — 2. The unit energy 
shell dF is given in terms of the function 



J s = Z(Jq) = — A 



i=i 1 



The volume of the region T is V r = n s/2 2 s-\ sr(s/2) UUi a i- 

The maximal volume V(Jn) of a cube touching the unit energy shell is given by 
V(J CT it) = IlUi a i (where J criM = jjfa). One thus finds the critical value 

_ 2^ S r( g /2) 

^ cnt > s ~~ n s/2 s s/2 y VL ) 
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Figure 4. Nodal count distributions for the harmonic oscillator (right column) and 
the cuboid (left column) for s — 2 (first row), s = 3 (second row), and s = 4 (third 
row). The black lines correspond to the limiting distributions -P s (£)- The red, green, 
and blue lines are numerically obtained histograms of the nodal count distribution at 
finite energy intervals: red line [E ,2E ], green line [2E ,AE ], blue line [4Eo,8Eq\. 
The chosen values for Eq are given in the corresponding graphs together with the 
chosen system parameters (the Hamiltonian) , and the bin size A£ that has been used 
for the histograms. In each case normalised nodal counts have been obtained for the 
lowest 50 to 80 million eigenfunctions. 

The insets magnify the graph near the critical cut-off value £ C rit- Overall the 
numerically obtained histograms are consistent with the (weak) convergence to the 
limiting distribution. 
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above which the limiting nodal domain distribution vanishes. Similarly to the harmonic 
oscillator the critical value and the limiting distribution do not depend on the detailed 
system parameters such as the side lengths. Indeed the limiting distribution may be 
written as 

™-^4*(-|f; i 

\ \ 1=1 / 1=1 j i=i 

For s = 2 this reduces to 

/ 2t2\-V2 o 

P 2 (fl=(l--i-J fOT ^<-- (63) 

See Figure 4 for graphs of the distribution for s = 2,3,4 together with numerically 
obtained histograms for finite energy intervals. 

5. Discussion 

We have derived an expression for the limiting nodal count distribution in the case 
where the wave equation is separable, and we have extracted some universal features of 
this distribution. While we have formally limited the scope with the assumptions (Al) 
to (A5) many standard examples as the harmonic oscillator or the particle in a cubic 
box obey all of these conditions. For some other examples which do not obey all of the 
conditions it is straight forward to generalize our derivations. For instance a particle in 
a spherical box does not obey assumption (A4) as action variables that correspond to 
angular momenta are not bounded from below. In this case the wave function separates 
in variables, some of which are cyclic. The derivation of a limiting distribution follows 
in full analogy once the expression for the nodal count in terms of quantum numbers is 
adapted and the Poisson summation is performed over the corresponding set of quantised 
action variables. Indeed most of our assumptions are purely technical and can be relaxed 
if necessary - though relaxing condition (A3) may imply that there are additional local 
maxima of the volume V(Jn) which may lead to further singularities within the support 
of the limiting nodal domain distribution. Also assumption (A5) that the Hamiltonian 
is a homogeneous function can be relaxed to a certain degree. Indeed it is only needed 
that the energy shell at high energies can be described asymptotically by a homogeneous 
function. 

It would certainly be interesting to compare our results to nodal count distributions 
of non-separable or wave-chaotic systems in dimensions larger than two. In any 
dimension one may try to obtain nodal counts numerically by using a corresponding 
adaptation of the Hoshen-Kopelman algorithm [16] and apply it to numerically obtained 
eigenfunctions. Berry's conjecture [17] states that highly excited chaotic eigenfunctions 
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can be simulated by a Gaussian random wave ensemble. In this ensemble the wave 
function is given by 



where nj are uniformly distributed on a unit (s — l)-sphere and the phases <f>j are 
equidistributed on [0, 2ir). On dimensional grounds one expects that the number of 
nodal domain in a given region of the random wave is proportional to the volume of the 
region. For two-dimensional random waves this has been checked and it is consistent 
with the critical percolation conjecture [18]. We tried to check this in three dimensions 
by finding the number of nodal domains of random waves inside a cube of side length 
a (at fixed wave number k — 1). The artificial boundary of the cube leads to nodal 
domains which intersect the boundary - indeed we have found that all nodal domains 
in our numerical approach were intersecting the boundary and that the nodal count 
is proportional to a 2 rather than a 3 . This scaling is expected on dimensional grounds 
for the number of nodal domains which intersect the boundary. However we were not 
able to increase the side length beyond a = 100 (about 16 wave lengths) on a standard 
desktop and we have just looked at a few hundred realisations. We cannot exclude that 
interior nodal domains (those which do not touch the artificial boundary) start to appear 
in much larger cubes and eventually dominate the nodal count. We can say however 
that any crossover from a 2 (boundary dominated) to a 3 (bulk dominated) would have 
to occur at considerably higher side lengths for which applying our numerical algorithm 
is beyond the power of standard desktop computers. 

Our numerical findings do confirm the basic expectation that the universality of a critical 
percolation model does not apply in the three dimensional case. For instance we find 
that the volume of the largest nodal domain scales linearly with the volume of the cube - 
a clear indication that one is inside the (non-universal) percolating regime (a non-trivial 
exponent is expected at the percolation transition). We hope that future research will 
shed more light on the nodal sets and nodal counts of wave-chaotic systems in dimensions 
s > 3 as well as in the corresponding random-wave models. 
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